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Mean field slider block models have provided an important entry point for understanding the 
behavior of discrete driven threshold systems. We present a method of constructing these models 
with an arbitrary frictional weakening function. This 'forced weakening' method unifies several 
existing approaches, and multiplies the range of possible weakening laws. Forced weakening also 
results in Abelian rupture propagation, so that an avalanche size depends only on the initial stress 
distribution. We demonstrate how this may be used to accurately predict the long-time event 
statistics of a simulation. 
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Complex spatial or temporal patterns often emerge 
when systems are driven from equilibrium and experi- 
ence instabilities. Coupled lattice maps have become a 
popular method of modelling such phenomena, seeking to 
capture the essential physics in a discrete form amenable 
to numerical simulation. In the study of earthquakes, 
simple slider-block models |^ ^ have long been employed 
to investigate the origin of magnitude-frequency scaling. 
Similar models are used to describe such diverse phe- 
nomena as pinned charge density waves flux lattices 
in type II superconductors , and creeping contact lines 

0- . 

Recently, analytic results for the avalanche size dis- 
tribution have been presented for mean-field slider-block 
models These models are characterized by homo- 

geneous, infinite range elastic interactions, and exhibit 
complex event histories and regimes of behavior. De- 
spite their simplicity, mean-field models remain sensitive 
to the choice of update rules and details of implemen- 
tation. This is especially evident in models that impose 
a particular form of simulated frictional weakening 
where different modes of behavior appear as the strength 
of the weakening is varied. 

Here we describe a method of constructing a mean- 
field threshold model where one may simulate arbitrary 
weakening laws under identical rules of evolution. This 
unifies the analysis of previously incompatible models, 
and provides more freedom in numerical simulation. This 
technique also results in Abelian rupture propagation, 
where the size of a simulated earthquake is uniquely de- 
termined from initial conditions, leading to a rigorous 
and implementation-independent analysis. 

The general slider-block model represents stick-slip 
motion along a fault plane with A'' ^ 1 discrete coor- 
dinates (or 'sites'). Each site i is assigned a slip deficit 
Ui, which measures the local distance from elastic equi- 
librium. The sites are pinned in place by frictional forces, 
and are subject to a restoring force (stress) proportional 
to their slip deficit. Internal disorder gives rise to an ad- 



ditional component of stress due to elastic interactions. 
The stress ai at a site i is related to the slip deficits 
through a linear constitutive relation 

ai = -KlUi - ^ Kij {ui - Uj) (1) 
j 

where and Kij are spring constants. If we impose uni- 
form (mean-field) interactions between all the elements, 
Kij — Kc/N, the above relation becomes 

ai = -KLUi - Kc{ui - (u)) (2) 

where (u) = ^i'^i will denote an average over all 
N sites in the model. We obtain a unitless expression by 
dividing by Kca, where a is a characteristic microscopic 
length. Defining the unitless slip deficit = u/a, stress 
a = 'aJ{Kca), and spring constant ratio — Kl/Kc, 
Eq. simplifies to 

a, = -{KR + l)q}, + {^). (3) 

For finite N we will refer to this as the near mean-field 
(NMF) model. Note that it is easy to invert Eq. (^ for 
the slip deficits in terms of stresses. 



Kn + 1 Kr{Kr + 1) ' ' 

so the configuration is uniquely determined by the pa- 
rameter Kfi and either the slip deficits or stresses alone. 

The model is slowly driven away from equilibrium 
by uniformly increasing the slip deficits. Eventually 
the stress at one site will surpass the maximum local 
frictional force and 'fail', sliding toward its equilibrium 
point. The motion of a failed site will change the mean 
slip deficit (</)), and produce a change in stress at other 
sites. If this change brings other sites to their threshold, 
they will also fail, producing an avalanche interpreted as 
a single event. 
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It is assumed that following the initiation of motion, 
the frictional force on a site will weaken, producing a 
transient dynamic instability. In discrete time we can- 
not model the dynamic slip or velocity of the site, but 
instead assign a residual stress cr^ at which the motion 
arrests. This is chosen from a probability distribu- 
tion independently for each failed site. Since slips occur 
instantaneously, we lose the interplay between a contin- 
uously evolving stress field and frictional force at a site. 
The behavior of dynamical models is known to strongly 
depend on the form of frictional weakening so the 
ability to include equivalent effects in discrete models 
would be advantageous. 

Since the NMF model is not dynamical we are only 
interested in large-scale features of its behavior that are 
independent of microscopic dynamics. Thus we are free 
to choose the simplest update rules that are consistent 
with the phenomenon of interest. In practice, we assume 
that a single site j reaches its stress threshold first. Since 
the physics will depend only on changes in stress, we may 
impose a uniform failure threshold by absorbing any 
variations into the residual stress distribution. The slip 
displacement — (jy'P — (j)"-'^ is related to the change in 
stress Ao-j = af - by Aj = -/\aj/{KR + 1- N''^). 
The motion of the site will change the mean slip deficit 
((/)) by Aj/N. We may view this as a transfer of stress 
from failing sites to all others. 

The above describes a series type dislocation where 
sites fail in sequence and the stress transfer occurs to 
other sites all at once. This makes it likely that any 
failing site (other than the single initiator) will have a 
stress slightly above the threshold, which subtly provides 
an order-of-failure dependence to the stress transfer. As 
a consequence, the exact stress transfer in simulation will 
depend on obscure factors like the order of iteration over 
sites. To eliminate this we must examine the stress trans- 
fer in more detail. 

Suppose that in the course of a (possibly ongoing) 
event there have been k block failures. Let {fc} repre- 
sent the set of indices of failed sites. Call k = k/N the 
fraction of failed sites. Then from Eq. (||) the change in 
stress for any stable site i is 

1 S 

= ((a/),, (5) 

where S = {Kr -I- 1 — 1/A^)~^, and is an average 
applied over failed sites. We call this the external stress 
transfer to signify that it applies to sites that are not 
part of the rupture. The quantity crj is the stress of site 
j at failure, which may be greater than . Note that 
the slip displacement Aj is only dependent on the stress 
drop at failure because pinning occurs immediately, and 
subsequent stress changes will not affect the slip. Since 
a-^ is typically very near and the aj^ are identically 
distributed random variables, the term in parenthesis is 



on average independent of k. Thus the external transfer 
grows linearly (on average) with the fraction of failed 
sites. 

In 'dynamic weakening' models |ic| ] the increased 
propensity for a failed site to slip further (due to a weak- 
ened pinning force) is simulated by imposing a lower 
threshold stress < for the duration of a single 
event. After a site fails it will receive stress transfer from 
subsequent failures, and thus may reach this lower thresh- 
old and fail again. Re-failing sites will contribute more 
to the external transfer and enhance the likelihood of 
continued rupture growth. 

This form of weakening is strange in that it first re- 
quires some failed sites to have their stress brought back 
up to the dynamical threshold, which will occur at some 
minimum rupture size. Following the onset of dynamical 
weakening the additional stress transfer typically results 
in a runaway event which fails every block in the sys- 
tem. This feature can be exploited to produce character- 
istic events which always occur once the minimum size is 
reached. 

One way to visualize the effects of weakening is to ex- 
amine the average stress of sites that have failed as a 
function of rupture size (Figure 1). Without weaken- 
ing, the average stress of failed sites (the average internal 
stress) will itself grow linearly like the external transfer 
(with half the slope). However, with dynamic weakening, 
all failed sites with stress > will fail again, putting a 
ceiling on the average internal stress. 

We claim that a more realistic approach would be to 
have failed sites shed a certain fraction of the stress they 
receive after failure. Implementing this 'fractional weak- 
ening' would involve re-computing the slip of all failed 
sites with each new failure. There is an easier way to ac- 
complish this if we note that the desired effect is to lower 
the slope of the average internal stress (AIS) function. 
There is a way to invert this relationship, such that the 
average internal stress function is given and the requisite 
slips and stress transfer computed as a result. In essence, 
in place of solving dynamical equations involving slip or 
velocity dependent friction, we can impose the effects of 
the weakening as they appear in a discrete time context. 
We call this approach 'forced weakening'. 

To perform this inversion we first observe the change 
in stress of a site j as it depends on fc, the number of 
failed sites. Let A*^ denote the slip displacement of site 
j after k failures. If it is nonzero it includes all block 
motion, including initial failure, additional failures from 
weakening, or continuous sliding. Consider the system of 
equations for the stress changes of the failed sites j G {k} 

= -(i^«-M)Aj- + l J2 (6) 
ieffc} 

The last line demonstrates the simple linear form of the 
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relationship between stress drops and slip displacements. 
This may be obtained via a matrix with diagonal ele- 
ments N^-^ — {Kii + 1) and off-diagonal elements N^^. 
This matrix may easily be inverted to obtain the slips A*^ 
in terms of the current stress drops Aa^ 

_ -{Kr + 1 - k/N)Aa^ ~ Afif 



-Act*-- 



(Kr + l)iKR + l- k/N) 
(k/N) {Aa'^)^ 



Kr + 1 {Kr + 1) [Kr + 1 - k/N) 



(7) 



This expression provides the slips A*^ that are necessary 
to generate a set of given stress drops. 

The slips are directly related to the external transfer, 
as in Eq. |. Summing over them yields a new expression 
for the external transfer in terms of the stress drops. 



Ac 



= 1 y A, = A /_A^'=\ 



je{fe} 
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Kr + \ {KR + \){KR + \-klN) 



{Kr + I^ k] 



(8) 



This is identical to Eq. (||) when fc = 1. However, the 
effective transfer coefficient S = S{k) — {Kr + 1 — k)~^ 
now grows with the rupture size. Making up for this is the 
fact that the stress drops are no longer frozen after failure 
and will now decrease as the rupture grows. Observe that 
to calculate the external transfer, we need not specify 
individual stress drops or slips, but we just need to know 
the average dynamic stress drop of all failed sites for a 
given rupture size. This evolving average stress drop is 
defined as 



k '^—^ 

3e{k} 



- - /(«) 



(9) 



where we have defined the average internal stress (AIS) 
function /(k) = (cr'^)^. The stress transfer now depends 
only on the initial stress and the AIS function, which 
are independent of failure order. Using this formulation, 
avalanches are Abelian. 

In numerical simulation, we must eventually assign an 
actual slip and/or residual stress to each failed site. Care 
must be taken to make results agree with the analogous 
forward simulation. For example, given a linear AIS func- 
tion /(k) ~ an + /3, we could assign the /c"* failed site a 
random residual stress (with mean /3) plus 2a{k — 1)/N. 
When using AIS functions with no forward equivalent, 
the method of assigning final stresses must be stated ex- 
plicitly. 

The forced weakening method has two main advan- 
tages. Practically, it allows the simulation of models with 
arbitrary weakening characteristics, some of which would 



not be obtainable with modified CA rules. Formally, the 
model is Abelian, so that the event size is a unique func- 
tion of the initial stress configuration. Using this fact we 
can seek an expression which will determine the event 
size given adequate information of the initial stresses. 

To accomplish this we first approximate the stress con- 
figuration with a continuous distribution Paix)^ where 
Pa-ix)d-X is the probability that a randomly selected site 
will have stress a between X; and x ~^ dx- There are 
two ways of obtaining a continuous probability distribu- 
tion from a discrete set of stresses: allow the number of 
sites N to become infinite, or to consider averages over 
a course grained time. In the latter case the continuous 
driving of the slip deficits will quickly fill in the spectrum 
of possible stresses. 

When considering a course grained time, we realize 
that the distribution is only an average, and that the ac- 
tual instantaneous distribution of stresses will fiuctuate. 
Similarly, the solution we obtain will be an average with 
statistical properties we would like to calculate. Under 
ordinary conditions, the distribution of stresses is sta- 



tistically stationary |11|, so the expected distribution of 
event sizes at any time will predict the long-time average 
behavior. 

To write down the next-event-size expression, it is con- 
venient to define the stress deficit and the 
cumulative distribution 



pMdx' 



(10) 



where ps(x) = P<y{'^^ ~ x)- Given this distribution the 
event size k is the solution to 



1 



^s(x)dx 



-cr^PjT^K + dK) = (11) 



where /(k) is again the AIS function. A full derivation 
and analysis of the results of this equation will be pub- 
lished elsewhere. For an example solution, suppose each 
stress in the system is chosen independently from a uni- 
form distribution between zero and . Then for large 
N , Eq. (|i^) determines that k is the first crossing of a 
1-d random walk constrained to return to zero at k = 1. 
Without that constraint, this produces a power law dis- 
tribution of event sizes with exponent —3/2. With the 
constraint in place, there is an upturn at large event sizes, 
which is observed in simulations (so long as Kr is low 
enough to allow such large events). A comparison of this 
analytical solution with a simulation is shown in Figure 
2. 

In summary, the forced weakening method extends 
the basic slider block model to include arbitrary weak- 
ening functions in an efficient discrete time simula- 
tion. Additionally, forced weakening results in Abelian 
avalanches which leads to a rigorous and implementation- 
independent analysis. We have presented a brief example 
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of this showing excellent agreement between theory and 
simulation. This approach should be equally applicable 
to related discrete threshold models. While the result- 
ing formalism is dependent on the mean-field character 
of the model, the numerical techniques may find wider 
use whenever an inversion like Eq. (0) is available. 
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FIG. 2. Comparison of a solution of Eq. (|1]J) for uniformly 
distributed stresses (a) with simulation results (b). The power 
laws have identical mean-field exponents of -3/2 and the pre- 
dicted upturn at large event sizes is evident in the simulated 
data. The parameters of the simulation are Kr = 0.01 and 
an average residual stress of 0.05. 




FIG. 1. Average internal stress (AIS) functions for several 
weakening laws. The AIS function measures the average stress 
of failed sites as a function of rupture size. 'Dynamic weaken- 
ing' as implemented in the literature places an upper bound 
on the AIS starting at the dynamic threshold . Fractional 
weakening, proposed here, sheds a fixed fraction of the inter- 
nal transfer to unfailed sites. In this figure we have taken the 
average residual stress to be 0.05. 
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